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Using the simple procedure, recently introduced, of dividing Gaussian matrices by 
a positive random variable, a family of random matrices is generated characterized 
by a behavior ruled by the generalized hyperbolic distribution. The spectral density 
evolves from the semi-circle law to a Gaussian-like behavior while concomitantly the 
local fluctuations show a transition from the Wigner-Dyson to the Poisson statistics. 
Long range statistics such as number variance exhibit large fluctuations typical of 
non-ergodic ensembles. 

I. INTRODUCTION 

In a recent paper [lj, the concept of disorder has been introduced in random matrix 
theory (RMT) to denote ensembles of random matrices generated by superimposing an 
external source of randomness to the fluctuations of the Gaussian ensembles. The method 
generates matrices which preserve unitary invariance though at the price of introducing 
correlations among matrix elements. Their statistical properties are amenable to ana- 
lytical derivation. Families of ensembles recently proposed in the literature fit into this 
category, for instance (3, ^ in the Wigner Gaussian case and, Refs. 0-0] in the case 
of Wishart matrices. Disordered ensembles can be considered as an instance of the so- 
called superstatistics[3] aside from the fact that their density distributions are normalized 
differently, as explained below. 

As a consequence of having operating an extra source of randomness, the ergodicity 
property of the standard ensembles, characterized by the equivalence between spectral 
averages taken along the spectrum and averages extracted from the spectra of a set 
of matrices is broken. This creates an ambiguity in measuring its collective spectral 
fluctuations although for those which concern individual eigenvalues such as extreme value 
statistics it remains well defined. This kind of statistics has been object of many studies^ 
in the last decades. In this direction, it was shown how the behavior of the largest 
eigenvalue described by the Tracy- Widom distribution^] is affected by disorder 0, 0]. 

Despite the great success the RMT Gaussian ensemble had and has (see [Io[ for recent 
reviews), it has features of an abstract mathematical model such as, for example, an 
eigenvalue density defined in a compact support. Parallel to its development, there 
was a search for models whose randomness would be closer to physical situations. This 
is the case of the two-body random ensemble (TBRE) proposed in the context of the 
shell model of nuclear physics [TO] and of the more general A;-body embedded Gaussian 
ensembles (EGE) proposed to many-body systems [121] . In them, the Hamiltonian model 
is made randomic by taking the strength of the residual interaction among particles from 
a Gaussian distribution. These ensembles show as main features: correlations among 
matrix elements, Gaussian eigenvalue density and non-ergodicityjH]. Nevertheless, their 



spectral fluctuations are supposed to follow, under restricted conditions that take into 
account the non-ergodicity ambiguity, the same fluctuations pattern of the Gaussian 
matrices, namely the Wigner-Dyson statistics [14]. Disordered ensembles may provide a 
simple way to understand features of these physical motivated ensembles. 

Families of disordered ensembles have, so far, been considered mostly in association 
with an external randomness governed by the gamma distribution 0, [sj] . This leads to 
ensembles in which statistics are dominated by power law behavior that implies strong 
non-ergodicity. It has also been considered the case of the inverse gamma distribution [5J] 
which contains singularities. Our purpose here is to study an ensemble whose auxiliary 
random variable is taken from a generalized inverse Gaussian (GIG) which acting on the 
Gaussian matrices generates disordered matrices distributed according to a generalized 
hyperbolic distribution (GH). This kind of distribution has been introduced in 197 7 by 
Barndorff-Nielsen 15( and has found since then many applications specially in finance {16]. 



II. DISORDERED ENSEMBLES 



We start by reviewing general aspects of the formalism of the disordered ensemble of 
Ref. |l[. The idea is to construct matrices by the relation 
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where Hy is a random matrix of dimension N with joint density probability distribution 
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and £ is a positive random variable with a normalized density probability distribution 
w(£) with average £ and variance cr|. In (|2J), f = N + f3N(N — l)/2 is the number of 
independent matrix elements, ft is the Dyson index /3 = 1,2,4 and V(x) is a confining 
potential which makes normalization 
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finite with respect to the measure dH = Yli=i^uY[j>iY[k=i'^^^ij- Rewritten as 
Hy = VsTC) Eq. ([I]) means that the joint distribution of the matrix elements of the 
disordered ensemble, is obtained from Eq. (J2J) by averaging over the £ variable, namely 
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Eq. (j3j) translates the relation ([T]) among the random quantities into a relation among 
their density distributions functions. Despite the resemblance of Eq. (jl]) with distri- 
butions constructed in the context of the superstatistics formalismj7], we remark that 
here the two distributions are integrated normalized while in the superstatistics case the 
normalization is performed after the integration. 



While Eq. flU) puts in evidence the correlations among matrix elements, the relation 
(PQ) which is expressed in terms of the random quantities themselves, gives a clearer 
picture of the model and, at the same time, provides an efficient algorithm to generate 
the disordered matrices. We also remark that in p], an alternative procedure has been 
described in which matrices of the ensemble are generated taking into account correlations 
among their elements (see next section). 

Turning now to eigenvalues and eigenvectors, we observe that by diagonalizing the 
matrices in flTJ, the relation 

is obtained, where D(£) and Dy are diagonal matrices whose elements are the respective 
eigenvalues of H(£) and Hy. Rewriting, as done with Eq. (CQ), Eq. (j5]) as Dy = 
D(OVUi, the relation 
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In (J7J), Kn is a normalization constant and J(x±,X2, ■■■Xn) is the eigenvalue part of the 
Jacobian of the transformation from matrix elements to eigenvalues and eigenvectors. 
From (jSJ), statistical measures of the generalized family can be calculated by weighting, 
with the distribution, the corresponding measures of the original ensemble with the 
eigenvalues multiplied by the factor 

Considering particular choices of the distribution w(£), the most used so far both, in 
the Gaussian and the Wishart cases, was the gamma distribution 

w{0 = ex P (-o^7r(0 (8) 

with variance = . Considering the case of the Wigner ensemble in which V(H) = 
H 2 , when (jSJ) is substituted into (BJ, the integrals are readily performed to give the joint 
distribution 

To illustrate the power-law behavior induced by the above gamma distribution, take out 
of the / independent elements, the density distribution of a given one denoted by h (see 
next section for the definition of the h variable), after integrating over the / — 1 others, 
we obtain 
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Since for large \h\ , pp(h] £) ~ 1/ |/i| 2 ^ +1 , the distribution (flQj) does not have moments of 
order superior to 2£. This fact makes the value £ = 1 critical since below it (fTOj) does 
not have a second moment j^j. We remark that, apart from the lack of independence, the 
marginal distribution of the matrix elements has the same kind of distribution, namely one 
with an asymptotic power-law behavior, as the ones of the ensemble of Levy matrices (l9|. 

Another distribution considered in |5| was the inverse distribution obtained by making 
£ = 1/v which leads to 

M,) = ^exp(-£)^- 2 . 
Again, the integral in <^ can be performed to give 
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which shows, as a drawback, the presence of a singularity at the origin. 

Our present purpose is to study the model generated by taking the disorder variable 
out of the generalized inverse Gaussian 
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where ^a(^) is the modified Bessel function of the third kind. This is a three parameter 
probability distribution where a and b are positive and A is real. It contains the above 
distributions (jSJ) and ffTTj) as special cases. In fact, in the limit a — > 0, apart from 
scaling, the gamma distribution, Eq. (jSj), is recovered while in the limit b — > the inverse 
distribution is obtained. The GIG first moment is 



and its variance is 



By taking the ratio 
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we have a better parameter to understand the effect of the disorder generated by the 
GIG distribution. First we remark that, by changing in ( Tl6l) the sign of A, the symmetry 
property K_ u (x) = K u (x) of the Bessel functions, makes the above ratio symmetric with 
respect to the value A = — 1. Secondly, for large orders, the Bessel functions can be 
approximated by the dominant term of their uniform approximation 17 
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Using this, the Bessel functions in (I16p can be approximated for large A as 
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where, in presence of one, the square of the term yab/X has been neglected. As a 
consequence, with a = b, the ratio ( TT6l) becomes practically constant when A increases as 
indeed can be seen in Fig. 1 . On the other hand, by fixing A and the product ab, the 
GIG also becomes more localized as the ratio b/a decreases, as illustrated in the same 
Fig. 1. From this analysis, we expect larger disorder for small values of the parameter 
A or large values of the ratio b/a, with the value A = — 1 playing a special role in this 
generalized family. 



III. HYPERBOLIC DISORDERED GAUSSIAN ENSEMBLE 

The Wigner ensemble is obtained by making V(H) = H 2 with H a Hermitian ma- 
trix. The number of independent matrix elements is of course f = N + (3N(N — l)/2 
and it is convenient to consider the set of / independent elements denoted by h n with 
n = 1,2,3,...,/. They are defined such that first iV ones are equal to the diagonal el- 
ements, that is, hi = Ha. The remain ones, the off-diagonal elements, are scaled as 
h n = \/2H^ with n — N + 1, N + 2, ...,/. In terms of these / independent variables, the 
joint distribution, Eq. (@J, together with Eq. ( Fl3l) becomes 
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Each h n variable is treated equally in ( 1201) and that is the reason to introduce this set 
of variables in which the inconvenient factor of \/2 multiplying the off-diagonal ones has 
been absorbed in the new variable. The integration over the extra variable £ makes 
the distribution in (120]) a correlated distribution of the set of variables, but the above 
distribution has the important property that removing by integration any one variable, 
the distribution of the remnant others has the same form, namely, we have 
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This property is called scale-invariant occupancy of phase space [18]. Of course, it is a 
manifestation of the uncorrelation among the matrix elements of the original ensemble. In 
this way, one can expect this property to hold in general for the disordered ensembles. On 



the other side, it is not expected to hold for generalizations of distributions constructed 
using the strict superstatistics recipe. 

As it occurred with the gamma distribution and its inverse, after substituting (fi~3~]) in 
dlj), the integrals are performed and we find the joint distribution of matrix elements and 
the distribution of a generic matrix element h to be given, respectively, by the generalized 
hyperbolic distributions 
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where with 1 < k < f, Xk = X + k/2. The scale- invariant property has been used to 
write Eq. (|23|) . Comparing these GH equations with the correspondent equations for the 
gamma disorder, we see that the GH ones contain in themselves the gamma ones but 
multiplied by the Bessel function that decays exponentially for large values of \h\. This 
leads to a more regular behavior as becomes clear from the expression 
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for the moment of arbitrary nth order of ( 1231) . From the symmetry property of the Bessel 
function, i.e. K_ u (x) = K v {x) we conclude that (I23p has all moments, resulting in milder 
fluctuations. 

Consider the identity 
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in which by definition each ratio gives the conditional probability p(hk\h\ 
the last element hk in the numerator argument once the preceding k — 1 ones have been 
determined. Using the property of scale-invariance, this probability can be written in 
terms of the Bessel functions as 
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A more compact expression is obtained doing the following 
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with 6 = b. With these definitions the probability distribution of one element once 
k — 1 other ones have already been sorted is 
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With k running from 2 to /, these equations form a set of univariate distributions which 
can sequentially be used to generate a matrix of the ensemble in which correlations among 
elements are taken into account. 

Integrating over all eigenvalues but one and multiplying by N, the eigenvalue den- 
sity is expressed in terms of the Wigner's semi-circle law density Pg{E) = y/2N — E 2 /-k 
of the Gaussian ensemble 3 as 
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and, since the ratio between the two averages is less than one, the disordered density is, 
at the origin, smaller than the semi-circle value. When \E\ — » oo, ^ max — > and therefore 
the interval of integration in f[3TJ|) collapses and a crude approximation to the integral is 
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which shows that the density has Gaussian tails. In Fig. 2, it is shown that with the 
parameter A fixed at the value —1, the density undergoes a transition from the semi-circle 
to a Gaussian-like shape which becomes more and more deformed as a = b decreases. 

By integrating the density, Eq. (1301) . from the origin to a value E, we obtain the 
function 
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where 
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The functions N(E) and N G (E) count the average number of eigenvalues in the interval 
(0, E) for the the disordered and the Gaussian ensembles, respectively. 

To measure spectral fluctuations, the dependence on the density is removed by the 
transformation x = N(E) which maps the actual spectrum into a new one with unit 
density. Starting with short range correlations, we define the spacing function that gives 
the probability E(0,s) that the interval (— | , | ) is empty. To calculate this function we 
integrate Eq. (J6]) over all eigenvalues outside the interval the ( — !,!) to obtain with 
s = 2N(l) 
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where for the Gaussian spacing we use the Wigner surmise for the GOE case (j3 = 1] 
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From this function the nearest neighbor distribution (NND) is derived as p(s) — d g (°' s ) 
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The above equations are exact and the spacings presented below in Fig. 3 were 
calculated with them. To understand these results we derive approximated spacings con- 
sidering that for large N, the two counting functions, N(E) and Nq(E) can be replaced, 
at the center of the spectrum, by N(E) = p(0)E and Nq(E) = p G (0)E. With these 
approximations, Eq. (1331) becomes 
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and its derivative F(0, s) = —dE(0, s)/ds 

f(o,s)= (&/a) !_ pw{0 l I .lic-n-v^^rne-*. (.39) 

2K x (Vab)p(0)^Jo 

We remark that using Eq. (|3T|) it can be verified that the condition F(0, 0) = 1 is 
satisfied. We can go further by expressing Eq. (|39|) in terms of Bessel functions as 
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Finally, in the same approximation the NND has the expression 
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Replacing now the Bessel function by its asymptotic behavior, assuming that as a function 
of s, its argument is large, we find the NND decays as 
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This decaying can present two limiting situations, first, if a and b increases, the second 
term inside the square root becomes smaller than one, in such a way that a Gaussian 
decay is obtained by expanding the square root. In the second situation, b decreases and 
makes the second term much greater than one such that it becomes the dominant term 
leading to an exponential decay. 

Therefore, as a function of its parameters, this model constitutes a family which locally 
describes intermediate cases between the Wigner-Dyson and the Poisson statistics. This 
is illustrated in Fig. 3 where the cumulative NND, F(s) = 1 — F(0, s), is plotted by the 
indicated values of the parameters. We see that as the density goes from the semi-circle 
to the Gaussian-like shape, concomitantly the spacing moves from the Wigner surmise to 
the Poisson distribution 1 — exp(— s). Further, at an intermediate point of the transition, 
which coincides with the density approaching a Gaussian shape, the cumulative NND 
approaches the so-called semi-Poisson distribution given by 1 — (1 + 2s) exp(— 2s)(20|. 

A statistics that measures spectral long range correlations and, also estimates non- 
ergodicity[l|, [2l[, is the variance S 2 (L) of the number of eigenvalues in the interval 
[-6/2,6/2] with L = 2N(6/2). In Q it is shown that this variance is given by 
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where Eg is Gaussian number variance. In the limit of large spectra, we again use the 
linear approximations of the two counting functions to obtain 
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where using the fact that Y!q(x) is a logarithmic smooth function of its argument, the 
integral was asymptotically performed. Eq. (j4"5|) shows that the effect of disorder in the 
number variance statistics is twofold: (i) it enhances it by rescaling the argument of the 
Gaussian expression and (ii) it introduces a superpoissonian quadratic term that affects 
large intervals. These effects are illustrated in Fig. 4, in which Eq. (1451) is calculated 



with parameters values which give a density and a spacing distribution closest to the 
Wigner-Dyson's ones in Fig. 2 and 3. 

To study the behavior of the largest eigenvalues in the limit of large matrix size N, 
one introduces the scaled variable 
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in terms of which the probability Eg,r(0, s) that the infinite interval (s, oo) is empty is 
known for the three symmetry classes [g]. For the unitary class, (3 = 2 
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where q(s) satisfies the Painleve II equation 
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q(s) ~ Ai(s) when s — > oo, 

where Ai(s) is the Airy function. For the orthogonal (/3 = 1) and symplectic (/3 = 4) 
classes the probabilities are given, respectively, by 
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The above equations give a complete description of the fluctuations of the eigenvalues 
at the edge of the spectra of the Gaussian ensembles. When perturbed, it has been 
shown [9j that for the three symmetry classes the above probabilities modify to 
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Equations (|53|) and f|54]) give a complete analytical description of the behavior of the 
largest eigenvalue once the function w(£) is chosen. In j^], asymptotic results, when the 
limit iV — > oo is taken, have been derived without specifying u>(£). To do this, it was 
assumed that the localization of w(£), given by the ratio cr w /£ was dependent on the 
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matrix size N. Considering the distribution w(£) independent of N keeping the ratio ; x 
fixed when the matrix size iV increases, for the three invariant ensembles, the function 
Eq,0 becomes a step function centered at £ = 2N£/t 2 . Therefore, in this regime, the 
probability distribution for the largest eigenvalue converges to 
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that shows a GIG distribution for the square of the largest eigenvalue. 
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IV. HYPERBOLIC DISORDERED WISHART MATRICES 



Taking now V = X^X with X being a rectangular matrix X of size (MxN), the 
Wishart ensemble is defined by the joint density distribution 



Pw(X) = 




where / = f3MN. From (|57|) . the elements of X are Gaussian distributed. The eigenvalues 
of the random matrix V = X^X have joint density distribution 
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This ensemble was introduced by the statistician J. Wishart [22| and plays an important 
role in statistical analysis j23| where V = X^X is a covariant matrix. More recently, it 
appeared associated in the chiral random ensemble 24| when the square of its eigenvalue 
is taken. 

It can be shown that, for large matrices, the density of these eigenvalues approaches 
the Marchenko-Pastur density [251] 
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in which, with c 
this density is 
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As in the Wigner case, disorder is introduced in the ensemble by defining new matrices 
through the relation 
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which replaced in (J57l) leads to an ensemble with joint density distribution of matrix 
elements 
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After substituting (fl3l) in (1621) . as it occurred in the Wigner case, integrals can be per- 
formed and we get 
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The joint distribution of eigenvalues is given by Eq. (JSJ) with P v replaced by Pyy. 
Eigenvalue measures are therefore obtained by averaging those of the Wishart matrices 
and, for instance, the eigenvalue density is given by 
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Of course, the above density extends beyond the Marchenko-Pastur limits x± as can be 
seen in Fig. 5. By taking the limit x — » oo, the integration interval collapses and a crude 
approximation to the integral gives the exponential decay 
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For the counting function we find after integrating the density 
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The Wishart ensemble belongs to a class of random ensembles whose spectral fluc- 
tuations are derived from a general formalism based on orthogonal polynomials, the 
Laguerre ones, and share the same universal statistiscal properties [261]. This universality 
in the case of the Wigner-Dyson statistics is manifested after mapping the eigenvalues 
into variables with density one. We consider two spacing functions which give the prob- 
ability of an interval to be empty: one at the bulk and the other at the inferior extreme 
of the spectrum. 



At the bulk, we take the interval (x — 8/2,x + 8/2) where x = M is the average position 
of the Marchenko-Pastur density. The spacing or gap function after introducing disorder 
is related to the unperturbed one by 
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where the spacing s is calculated with Eq. (|67|) as 
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Comparing f )68|) and fl35|) . we conclude that these spacings have the same behavior. 

To study the behavior of the smallest eigenvalues in the limit of large matrix size N, 
one introduces the scaled variable s = x/4N in terms of which the probability E W b (0, s) 
that the interval (0, s) is empty has been derived in Refs. 27, 28|| for (3 = 2 and 29] for 
(3 = 1. When perturbed, the above probabilities modify to 



£7(0, 



(b/a 



I A/2 



2K x (Vab) 



<& A_1 exp 




(70) 



which shows how the smallest eigenvalue distribution is deformed in the presence of a 
hyperbolic disorder. This probability becomes simple when M = N with the Marchenko- 
Pastur density diverging at the origin. In this case the probability E w ^ (0, s) assume 
simple exponential form particularly for (3 = 1. 



V. CONCLUSION 

We have investigated the effect of superimposing an extra source of randomness gov- 
erned by a generalized inverse Gaussian to the Gaussian fluctuations of the Wigner and 
the Wishart ensembles. The result is an ensemble of random matrices ruled by the gen- 
eralized hyperbolic distribution which contains as particular cases disordered ensembles 
previously studied. The spectral density and short range statistics like spacing distribu- 
tion show a transition from Wigner-Dyson to Poisson statistics, approaching universal 
critical statistics, namely semi-Poisson statistics, at an intermediate point of the transi- 
tion. However, differently from the short range, long range statistics of the hyperbolic 
ensemble show a superpoissonian behavior with large fluctuations. The combination of 
semi-Poisson at short range and super-Poisson at long range have been observed in the 



non-ergodic embedded Gaussian models of many-body systems [131]. Therefore the present 
hyperbolic model is more suited to give a simple way to model features associated with 
the non-ergodicity of physical motivated ensembles. 

This work is supported by a CAPES-COFECUB project and by the Brazilian agencies 
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Figure Captions 

Fig. 1 The GIG distribution, Eq. ( IT3|) . is plotted for the indicated values of the 



parameters: in a) the parameters a and b are kept constants while in b) A = 5 and the 
ratio a/b is varied keeping ab = 4. 

Fig. 2 The transition of the density from the semi-circle to a Gaussian-like for the 
indicated values of the parameters. 

Fig. 3 The transition of the spacing F(s) from Wigner-Dyson to Poisson statistics is 
shown for the same values of the parameters of Fig. 2. 

Fig. 4 The nearest-neighbor distribution (NND) calculated with the approximated 
equation (1421) for the same values of the parameters of Figs. 2 and 3. 

Fig. 5 The number variance with small level of disorder is shown and compared with 
the GOE, the Poisson and the semi-Poisson predictions. 

Fig. 6 The density of the disordered Wishart matrices of size M=80 and N=40 is 
compared with the Marchenko-Pastur density for the indicated values of the parameters. 



